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We investigate the stability of possible order parameter configurations in clean layered heterostruc- 
tures of the SFS...FS type, where S is a superconductor and F a ferromagnet. We find that for 
most reasonable values of the geometric parameters (layer thicknesses and number) and of the ma- 
terial parameters (such as magnetic polarization, wavevector mismatch, and oxide barrier strength) 
several solutions of the self consistent microscopic equations can coexist, which differ in the arrange- 
ment of the sequence of "0" and "7r" junction types (that is, with either same or opposite sign of 
the pair potential in adjacent S layers). The number of such coexisting self consistent solutions 
increases with the number of layers. Studying the relative stability of these configurations requires 
an accurate computation of the small difference in the condensation free energies of these inhomo- 
geneous systems. We perform these calculations, starting with numerical self consistent solutions of 
the Bogoliubov-de Gennes equations. We present extensive results for the condensation free energies 
of the different possible configurations, obtained by using efficient and accurate numerical methods, 
and discuss their relative stabilities. Results for the experimentally measurable density of states are 
also given for different configurations and clear differences in the spectra are revealed. Comprehen- 
sive and systematic results as a function of the relevant parameters for systems consisting of three 
and seven layers (one or three junctions) are given, and the generalization to larger number of layers 
is discussed. 

PACS numbers: 74.50+r, 74.25.Fy, 74.80.Fp 



I. INTRODUCTION 

A remarkable manifestation of the macroscopic quan- 
tum nature of superconductivity is seen in the description 
of the superconducting state by a complex order param- 
eter with an associated phase, <p, which is a macroscopic 
quantum variable. For composite materials comprised 
of multiple superconductor (S) layers separated by non- 
superconducting materials, the phase difference A0 be- 
tween adjacent S layers becomes a very relevant quantity. 
For the case where a nonmagnetic normal metal is sand- 
wiched between two superconductors, it is straightfor- 
ward to see that the minimum free energy configuration 
corresponds to that having a zero phase difference be- 
tween the S regions, in the absence of current. The situ- 
ation becomes substantially modified for superconductor- 
ferromagnet-superconductor (SFS) junctions, where the 
presence of the magnetic (F) layers leads to spin-split 
Andree-\*i states and to a spatially modulated 2 * 3 - order 
parameter that can yield a phase difference of A</> = tt 
between S layers. These are the so-called tt junctions. 
Junctions of this type can occur also in more compli- 
cated layered heterostructures of the SFSFSF.. type, 
where the relative sign of the pair potential A(r) can 
change between adjacent S layers. 

Continual improvements in well controlled deposi- 
tion and fabrication techniques have helped increase the 
experimental implementations of systems containing tt 
junctionsi 4 ' 5 ' 6 ' 7 : 8 ' 9 ' 10 : 11 ! 12 ^ 3 Possible applications to de- 
vices and to quantum computing, as well as purely sci- 



entific interest, have stimulated further interest in these 
devices. The pursuit of the tt state has consequently 
generated ample supporting data for its existence and 
properties: The observed nonmonotonic behavior in the 
critical temperature was found to be consistent with the 
existence of a tt stated The domain structure in F is ex- 
pected to modify the critical temperature behavior how- 
ever, depending on the applied field* 7 - The ground state 
of SFS junctions has been recently measured^ and it 
was found that or tt coupling existed, depending on the 
width cIf of the F layer, in agreement with theoretical 
expectations. Similarly, damped oscillations in the criti- 
cal current Ic as a function of cIf suggested also a to 
7r transition^ The reported signature in the character- 
istic Ic curves also indicated a crossover from the to 
tt phase in going from higher to lower temperatures i& 
Furthermore, the current phase relation was measured, 12 
demonstrating a re-entrant Ic with temperature varia- 
tion. 

A good understanding of the mechanism and robust- 
ness of the 7r state in general is imperative for the further 
scientific and practical development of this area. The 
tt state in SFS structures was first investigated long 
ago^* In general, the exchange field in the ferromag- 
net shifts the different spin bands occupied by the corre- 
sponding particle and hole quasiparticles. This splitting 
determines the spatial periodicity of the pair amplitude 
in the F layeoi^ and can therefore induce a crossover 
from the state to the tt state as the exchange field ho 
varies, 16 - or as a function of (If- The Josephson critical 
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current was found to be nonzero at the to ir transition, 
as a result from higher harmonics in the current-phase 
relationship. 1 '' It was found that a coexistence of sta- 
ble and metastable states may arise in Josephson junc- 
tions, which was also attributed to the existence of higher 
harmonics^ If the magnetization orientation is varied, 
the 7r state may disappear, 19 in conjunction with the 
appearance of a triplet component to the order param- 
eter. A crossover between and ir states by varying 
the temperature was explained within the context of a 
Andreev bound state model that reproduced experimen- 
tal findings^ In the ballistic limit a transition occurs if 
the parameters of the junction are close to the crossover 
at zero temperature^ An investigation into the ground 
states of long (but finite) Josephson junctions revealed 
a critical geometrical length scale which separates half- 
integer and zero flux states^ The lengths of the individ- 
ual junctions were also found to have important impli- 
cations, as a phase modulated state can occur through a 
second order phase transition^ 

Nevertheless, little work has been done that studies 
the stability of the n state for an SFS junction from 
a complete and systematic standpoint, as a function of 
the relevant parameters. This is a fortiori true for lay- 
ered systems involving more than one junction, each of 
which can in principle exist in the or the 7r state, or 
for superlattices. There are several reasons for this. The 
existence or absence of n junction states is intimately 
connected to the spatial behavior of A(r) and of the pair 
amplitude F(r) (which oscillates in the magnetic layers), 
and thus the precise form of these quantities must be cal- 
culated self-consistently, so that the resulting F(r) cor- 
responds to a minimum in the free energy. An assumed, 
non self consistent form for the pair potential, typically 
a piecewise constant in the S layers, often deviates very 
substantially from the correct self consistent result, and 
therefore may often lead to spurious conclusions. Indeed, 
the assumed form may in effect force the form of the fi- 
nal result, thereby clearly leading to misleading results. 
Self-consistent approaches, despite their obvious superi- 
ority, are too infrequently found in the literature primar- 
ily because of the computational expense inherent in the 
variational or iterative methods necessary to achieve a 
solution. A second problem is that investigation the rela- 
tive stability of self-consistent states requires an accurate 
calculation of their respective condensation energies, and 
that, too, is not an easy problem. 

In this paper, we approach in a fully self consistent 
manner the question of the stability of states containing 
7r junctions in SFSF...S multilayer structures. As has al- 
ready been shown in trilayers^i it can be the case that, 
for a given set of geometrical and material parameters, 
more than one self consistent solution exists, each with 
a particular spatial profile for A(r) involving a junction 
of either the or tt type. It will be demonstrated be- 
low that this is in fact a very common situation in these 
heterostructures: one can typically find several solutions, 
all with a negative condensation energy, i.e., they are all 



stable with respect to the normal state. Thus a careful 
analysis is needed to determine whether each state is a 
global or local minimum of the condensation free energy. 
This determination is, as we shall see, very difficult to 
make from numerical self-consistent results, because it 
requires a very accurate computation of the free energies 
of the possible superconducting states, from which the 
normal state counterpart must then be subtracted. This 
subtraction of large quantities to obtain a much smaller 
one makes the problem numerically even more challeng- 
ing than that of achieving self-consistency, since although 
the number of terms involved is the same, the numerical 
accuracy required is much greater. Until now, this has 
been seen as a prohibitive numerical obstacle. Removing 
this obstacle involves a careful analysis and computation 
of the eigenstates for each state configuration, that is, 
the energy spectrum of the whole system. 

The numerical method we will discuss and implement 
overcomes these difficulties, and therefore enables us to 
determine the relative stability of the different states in- 
volved, as we shall see, for a variety of F/S multilayer 
structure types, and broad range of parameters. The 
condensation energies for the several states found in a 
fully self consistent manner, are accurately computed as a 
function of the relevant parameters. The material param- 
eters we investigate include, interfacial scattering, Fermi 
wave vector mismatch, and magnetic exchange energy, 
while the geometrical parameters are the superconduc- 
tor and ferromagnet thicknesses, and total number of 
layers. We will see that as the number of S layers in- 
creases, the number of possible stable junction configu- 
rations correspondingly increases. Our emphasis is on 
system sizes with S layers of order of the superconduct- 
ing coherence length £o , separated by nanoscale magnetic 
layers. In order to retain useful information that de- 
pends on details at the atomic length scale, it is neces- 
sary to go beyond the various quasiclassical approaches 
and use a microscopic set of equations that does not av- 
erage over spatial variations of the order of the Fermi 
wavelength. This is particularly significant in our mul- 
tiple layer geometry, where interfering trajectories can 
give important contributions to the quasiparticle spectra, 
owing to the specular reflections at the boundaries, and 
normal and Andreev reflections at the interfaces^ The 
influence of these microscopic phenomena is neglected in 
alternative approaches involving averaging over the mo- 
mentum space governing the quasiparticle paths. Thus, 
our starting point is fully microscopic the Bogoliubov de- 
Gennes (BdG) equations,? 6 which are a convenient and 
physically insightful set of equations that govern inhomo- 
geneous superconducting systems. It is also appropriate 
for the relatively small heterostructures we are interested 
in, to consider the clean limit. 

The outline of the paper is as follows. In Sec. [D] we 
write down the relevant form of the real-space BdG equa- 
tions, and establish notation. After introducing an ap- 
propriate standing wave basis, we develop expressions for 
the matrix elements needed in the numerical calculations 
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of the quasiparticle amplitudes and spectra. The itera- 
tive algorithm which embodies the self-consistency proce- 
dure is reviewed. We then explain how to use the self con- 
sistent pair amplitudes and quasiparticle spectra to cal- 
culate the free energy, as necessary to distinguish among 
the possible stable and metastable states. In Sec. IIIII we 
illustrate for several junction geometries the spatial de- 
pendence of the pair amplitude F(r), which is a direct 
measure of the proximity effect and gives a physical un- 
derstanding of the various self consistent states we find. 
Wc display this quantity as a function of the important 
physical parameters. The stability of systems containing 
various numbers of it junctions is then clarified through 
a series of condensation energy calculations that again 
take into consideration the material and geometrical pa- 
rameters mentioned above. Of importance experimen- 
tally is the local density of states (DOS), which we also 
illustrate for certain multilayer configurations. We find 
differing signatures for the possible configurations that 
should make them discernible in tunneling spectroscopy 
experiments. To conclude, in Sec. II VI we summarize the 
results. 



II. METHOD 

A. Basic equations 

In this section we briefly review the form that the 
Bogoliubov-de Gennes (BdG) equations^ take for the 
S/F multilayered heterostructures we study. Additional 
details can be found in Refs. I2I27II2I The BdG equa- 
tions are particularly appropriate for the investigation of 
the stability of layered configurations in which the pair 
amplitude may or may not change sign between adjacent 
superconducting layers. These are conventionally called 
"0" or "tt" junction configurations respectively. 

We consider three-dimensional slab-like heterostruc- 
tures translationally invariant in the x — y plane, with 
all spatial variations occurring in the z direction. The 
heterostructure consists of superconducting, S, and fer- 
romagnetic, F, layers. Examples are depicted in Fig. Q 
The corresponding coupled equations for the spin-up and 
spin-down quasiparticle amplitudes (u^v^) then read 



1 d 2 

~2m~~dz 2+£i - ~ Ep ^ + U W ~ h o( z )Wn( z ) + A ( z ) v k( z ) = e»<0) 
1 d 2 

"0 — TP! + 6 ^- E ^ z ) + U ( z ) + h °( z )K( z ) + A ( Z ) U U Z ) = e n vi(z) : 
Am oz z 



(la) 
(lb) 



where e± is the kinetic energy term corresponding to 
quasiparticles with momenta transverse to the z direc- 
tion, e n are the energy eigenvalues, A(z) is the pair po- 
tential, and U(z) is the potential that accounts for scat- 
tering at each F/S interface. An additional set of equa- 
tions for and v\ can be readily written down from 
symmetry arguments, and thus is suppressed here for 
brevity. The form of the ferromagnetic exchange en- 
ergy ho(z) is given by the Stoner model, and therefore 
takes the constant value ho in the F layers, and zero else- 
where. Other relevant material parameters are taken into 
account through the variable bandwidth Ep(z). This 
is taken to be E F (z) = E FS in the S layers, while 
in the F layers one has Ep(z) — Efm so that in 
these regions the up and down bandwidths are respec- 
tively Ef\ = Efm + ho, and Ef\, — Efm — ho. The 
dimensionless parameter /, defined as / = ho /Efm, 



conveniently characterizes the magnets' strength, with 
1 = 1 corresponding to the half metallic limit. The ratio 
A = Efm/Efs = (&fm /kFs) 2 describes the mismatch 
between Fermi wavevectors on the F and S sides, assum- 
ing parabolic bands with kps denoting the Fermi wave 
vector in the S regions. 

The spin-splitting effects of the exchange field coupled 
with the pairing interaction in the 5* regions, results in a 
nontrivial spatial dependence of the pair potential, which 
is further compounded by the normal and Andreev scat- 
tering events that occur at the multiple S/F interfaces. 
When these complexities are taken into account, one gen- 
erally cannot assume any explicit form for A(z) a priori. 
Thus, when solving Eqs. QJ, the pair potential must be 
calculated in a self consistent manner by an appropriate 
sum over states: 



A(*) = E / d£ ^- R( z >n( z ) + tanh( £ „/2T), (2) 

I 

where N(0) is the DOS per spin of the superconductor in the normal state, d is the total system size in the z 
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(b) 
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FIG. 1: Examples of the two types of multilayer geometries 
for the heterostructures examined in this paper. The system 
has a total thickness d in the z direction, and the F layers have 
thickness cLf. The general patterns shown hold for structures 
with an arbitrary odd value of the number of layers, Nl- The 
seven layer case is displayed. In pattern (a) the thicknesses 
of each S layer is ds, while in (b) the two outer S layers have 
thickness ds, and the inner ones have thickness 2ds (see text). 



direction, T is the temperature, lod is the cutoff "Debye" 
energy of the pairing interaction, and g(z) is the effective 
coupling, which we take to be a constant g within the 
superconductor regions and zero elsewhere. 

The presence of interfacial scattering is expected to 
modify the proximity effect. We assume that every S/F 
interface induces the same scattering potential, which we 
take of a delta function form: 



where zi is the location of the interface and H is the 
scattering parameter. It is convenient to use the dimen- 
sionless parameter Hb = mH/kps to characterize the 
interfacial scattering strength. 

An appropriate choice of basis allows Eqs. (JIJ to be 
transformed into a finite 2N x 2N dimensional matrix 
eigenvalue problem in wave vector space: 



H+ D 
D H~ 

T 



,1 



zN' u nl> ■ 



iN 



(4) 



), are the ex- 



where 'F^ = (u 
pansion coefficients associated with the set of orthonor- 
mal basis vectors, u\{z) 



;(*) 



2 / d E^=i4<2 sin(fc 9 z) 



and 

The longitudinal mo- 



U{z h z)=HS(z- Zl ) 



(3) 



mentum index k q is quantized via k q = q/ird, where q 
is a positive integer. The label n encompasses the index 
q and the value of ex- The finite range of the pairing 
interaction ojd, implies that N is finite. In our layered 
geometry submatrices corresponding to different values 
of e± are decoupled from each other, so one considers 
matrices labeled by the q index, for each relevant value 
of e±_ . The matrix elements in Eq. Q depend in general 
on the geometry under consideration, and are given for 
two specific cases in the subsections below. 



B. Identical superconducting layers 

The first type of structure we consider is one consisting 
of alternating S and F layers, each of width ds and dp 
respectively. This geometry is shown in Fig. ^a) for the 
particular case of Nl = 7. For a given total number of 
layers (superconducting plus magnetic) Nl, we have in 
this case for the interfacial scattering: 

(N L -l)/2 

U{z)= [U{i{ds+d F ),z)+U{{i{ds+d F )-d F )),z] 

i=l 

where U(zi,z) is given in Eq. (|3J). The matrix elements 
Hq q , and H~ q , in Eq. (0} are compactly written for this 
geometry as 



Hi = 

2m 



Efs 



2H 
d 1 2 



N, 



- A(2q) 



-E 



dp 
~~d 



Nr 



B(2q) 



- B(2q) 



2H 



K' = -^-[Ml - </) - A (Q + <?')] + I^FT - E FS ][B{q - q') - B(q + q' 



(6a) 
(6b) 



5 



where, 



A(q) = cos 



d F irq\ 



(N L -l)/2 



2d 



B(q) 



=1 

— > cos 



2 sin 



7rg 



7rg 

7rq 
2d 



(2d s i + d F (2i - 1)) 
(d s (2i-l)+2d F (i-l)) 



The matrix elements H , are similarly expressed in term of the coefficients A(q) and B(q), 



— 

H„„ = —— 6 

qq 2m 



Efs 



2H 
~d 

ds (N L + \ 
d 



N L ~1 
2 

- B(2q) 



A(2q) 



E 



dp (N L -1 
d 



B{2q) 



2H, 



H~, = —[A(q - q') - A(q + q 1 )] - [E Fl - E FS ][B(q - q') - B(q + q% 



r 



q^q'- 



(7a) 
(7b) 



(8a) 
(8b) 



The D qq i in the off-diagonal part of the left side of 
Eq. J2| arise from an integral over A(z), which scatters a 
quasiparticle of a given spin into a quasihole of opposite 
spin. One has: 



D qql = - 



dz sin(fc g z) A(z) sin(k'z). 



(9) 



It is straightforward to write also2i2& the self consistency 
equation in terms of matrix elements. 



C. Half-width superconducting outer layers 



the same width ds- We are also interested in investi- 
gating structures where the inner 5* layers are twice as 
thick (2ds ) as the outer ones (see Fig. ^b)) while the 
F layers remain all of the same width. This case is of 
interest because, roughly speaking, the inner layers, be- 
ing between ferromagnets, should experience about twice 
the pair-breaking effects of the exchange field than do the 
outer ones. Therefore, the results might depend on 
more systematically, particularly for relatively small Nl, 
if the witdth of the outer layers is halved. This has been 
found to be the case in some studies^ of the transition 
temperature in thin layered systems. 



The previous subsection outlined the details needed to 
arrive at the matrix elements when the S layers are of 
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A slight modification to the previous results yields the 
following form of the scattering potential U(z), 



{Nl-\)/2 

U{z)= [U{(i(2d s + d F )-ds),z) + U((i(2d s + d F )-ds-d F ),z)] 



(10) 



The matrix elements H^ q , are now expressed as 
2H , 



Jit, : ~ [A(q -q , )-A(q + q')} + [E F] - E FS ] [B(q + q>) - B(q - q')} , q^q' 
~N L - 1 



TT+ - ^ -x- c x 2H 
Hqq -2^ +£± + — 



A(2q) 



E 



FT 



dp (N L - 1 
d 



B(2q) 



E F s 



^(N L -l)+B(2q) 



where the coefficients A{q) and B(q) now read, 

(N L -l)/2 



A(q) = cos 



d F nq 
2d 



E 



cos 



^(2z - l)(2d s + d F ) 



B(,) = ^Mtanfe^ 



Tvq 



\ 2d J 



(11a) 



(lib) 



(12) 
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In a similar manner, The matrix elements H , are written as 



H qq = 



k 2 
K q 

2m 



■ e± 



2H 



N L -1 



A(2q) 



dp ( N L -1 
d V 2 



B(2q) 



E 



FS 



2H, 



(N L -l) + B(2q) 



= —[A(q - q') - A(q + </)] - [E Fl - E FS ][B(q + q') - B(q - q') 



q + q' 



(13a) 



(13b) 



The matrix elements of D are as in the previous sub- 
section, and the self consistent equation can be similarly 
rewritten. 



D. Spectroscopy 

Experimentally accessible information regarding the 
quasiparticle spectra is contained in the local density of 
one particle excitations in the system. The local density 
of states (LDOS) for each spin orientation is given by 

n 

(14) 

where a =|, J, and /'(e) = df/de is the derivative of the 
Fermi function. 

As discussed in the Introduction, the condensation free 
energies of the different self consistent solutions found 
must be compared^* in order to find the most stable 
configuration, as opposed to those that are metastable. 
While for homogeneous systems this quantity is found in 
standard textbooks (22*21. the case of an inhomogeneous 
system is more complicated. We will use the convenient 
expression found in Ref. [3^ for the free energy T: 

^=-2T£ln[2cosh(^)]+^z^, (15) 

where the sum can be taken over states of energy less 
than lod- For a uniform system the above expression 
properly reduces to the standard textbook result^ The 
corresponding condensation free energy (or, at T = 0, 
the condensation energy) is obtained by subtracting the 
corresponding normal state quantity, as discussed below. 
Thus, in principle, only the results for A(z) and the ex- 
citation spectra are needed to calculate the free energy. 
As pointed out in Ref. |3 a numerical computation of 
the condensation energies that is accurate enough to al- 
low comparison between states of different types requires 
great care and accuracy. Details will be given in the next 
section. 



III. RESULTS 

As explained in the Introduction, the chief objective of 
this work is to study the relative stability of the different 
states that are obtained through self-consistent solution 
of the BdG equations for this geometry. These solutions 
differ in the nature of the junctions. Each junction be- 
tween two consecutive S layers can be of the "0" type 
(with the order parameter in both S layers having the 
same phase) or of the 'V type (opposite phase) . As the 
number of layers, and junctions, increases, the number of 
order parameter (or junction) configurations which are 
in principle possible increases also. As we shall see, for 
any set of parameter values (geometrical and material) 
not all of the possible configurations are realized: some 
do not correspond to free energy local minima. Among 
those that do, the one (except for accidental degenera- 
cies) which is the absolute stable minimum must be deter- 
mined, the other ones being metastable. We will discuss 
these stability questions as a function of the material and 
geometrical properties, as represented by dimensionless 
parameters as we shall now discuss. 



A. General considerations 

Three material parameters are found to be very im- 
portant: one is obviously the magnet strength /. We will 
vary this parameter in the range from zero to one, that 
is, from nonmagnetic to half-metallic. The second is the 
wavevector mismatch characterized by A = (Ufm /&fs) 2 - 
The importance of this parameter can be understood by 
considering that, even in the non self consistent limit, 
the different amplitudes for ordinary and Andreev scat- 
tering depend strongly on the wavevectors involved, as 
it follows from elementary considerations. We will vary 
A in the range from unity (no mismatch) down to 1/10. 
We have not considered values larger than unity as these 
are in practice infrequent. The third important dimen- 
sionless parameter is the barrier height Hb defined below 
Eq. J2J). This we will vary from zero to unity, at which 
value the S layers become, as we shall see, close to being 
decoupled. We will keep the superconducting correlation 
length fixed at kps£,o = So = 100. This quantity sets 
the length scale for the superconductivity and therefore 



7 



can be kept fixed, recalling only that, to study the ds 
dependence, one needs to consider the value of ds/t;o- 
Finally the cutoff frequency uir> can be kept fixed (we set 
u)d = 0.04i?Fs) since it sets the overall energy scale and 
we are interested in relative shifts. The dimensionlcss 
coupling constant gN(0) can be derived from these quan- 
tities using standard relations. In this study we will focus 
on very low temperatures limit, fixing T to T = O.OlXjS 
(where Tq is the bulk transition value) . The geometrical 
parameters are obviously the number of layers Nl, and 
the thicknesses dp and d$- We will consider two exam- 
ples of the first: Nl = 3 and Nl — 7. For the larger 
value we will study both of the geometries in panels (a) 
and (b) of Fig.Q The thicknesses, usually expressed in 
terms of the dimensionlcss quantities Ds = kpsds and 
Dp = kpsdp 1 will be varied over rather extended ranges. 

As we study the effect of each one of these parameters 
by varying it in the appropriate range, we will be hold- 
ing the others constant at a certain value. Unless other- 
wise indicated, the values of the parameters held constant 
will take the "default" values I = 0.2,L» S = 100 = S , 
Dp = 10, A = 1, and Hb = 0. One important derived 
length is (fej — fcjj , where k-\ and k± are the Fermi 
wavevectors of the up and down magnetic bands. As 
is well knowniSiSI, this quantity determines the approx- 
imate spatial oscillations of the pair amplitude in the 
magnet. In terms of the quantities I and A we can de- 
fine: 



3 2 = k F s{k^-ki) 1 = 



1 



1 



(16) 



At I = 0.2 and A = 1 one has S 2 = 4.97 increasing 
to 15.7 at A = 0.1. This motivates our default choice 
Dp = 10. 

The numerical algorithm used in our self consistent cal- 
culations follows closely that of previous developed codes 
used in simpler geometries £12^ There are however some 
extra complexities that arise for the larger multilayered 
structures studied here, and from the increased number 
of self consistent states to be analyzed. As usual, one 
must assume an initial particular form for the pair po- 
tential, to start the iteration process. This permits a 
straightforward diagonalization of the matrix given in 
Eq. (0J for a given set of geometrical and material pa- 
rameters, for each value of the transverse energy e±. The 
initial guess of A(z) is always chosen as a piecewise con- 
stant ±Ao, where Ao is the zero temperature bulk gap, 
and the signs depend on the possible configuration be- 
ing investigated (see below) . Self consistency is deemed 
to have been achieved when the difference between two 
successive A(z)'s is less than 10 _5 Ao at every value of 
z. The minimum number of e± variables needed for self 
consistency is around N± = 500 different values of e±. In 
practice however, use of a value close to this minimum is 
insufficient to produce smooth results for the local DOS. 
Therefore, we first calculate A(z) self consistently us- 
ing N± = 500, after which iteration is continued with 
N±_ increased by a factor of ten. The computed spectra 



then summed according to Eq. (|14|) are smooth and fur- 
ther increases in N± produce no discernible change in the 
outcome. This two-step procedure leads to considerable 
savings in computer time. The convergence properties 
and net computational expense obviously depend also on 
the dimension of the matrix to be diagonalized, dictated 
by the number of basis functions, N, which scales lin- 
early with system size kpsd. Thus, the computational 
time then increases approximately as (kp$d) 2 , which can 
be a formidable issue. For the largest structures consid- 
ered in this paper, the resultant matrix dimension is then 
2N x 2N ~ 2000 x 2000 for each e x . 

The number of iterations needed for self consistency 
depends on the initial guess. If the self consistent state 
turns out to be composed of junctions of the same types 
as the initial guess, as specified by the signs in the ±Ao, 
then forty or fifty iterations are usually sufficient. But a 
crucial point is that, as we shall see, not all of the initial 
junction configurations lead to self consistent solutions 
of the same type. Since the self-consistency condition is 
derived from a free energy minimization condition, this 
means that some of the initial guesses do not lead to 
minima in the free energy with the same junction config- 
uration as the initial guess, and thus that locally stable 
states of that type do not exist, for the particular pa- 
rameter values studied. In this case, the number of iter- 
ations required to get a self consistent solution increases 
dramatically, since the pair potential has to reconfigure 
itself into a much different profile. The number of itera- 
tions in those cases can exceed 400. 

The computation of the condensation free energies of 
the self consistent states found is still a difficult task, even 
after the spectra are computed: considering for illustra- 
tion the T = limit, recall^ that the bulk condensation 
energy, given by 



E b = -(1/2)A(0)A^ 



(17) 



represents a fraction of the energy of the relevant elec- 
trons of order (A/wb) 2 , a quantity of order 10~ 2 or 
less. The condensation free energies of the inhomoge- 
neous states under discussion are usually much smaller 
(by about an order of magnitude as we shall see) than 
that of the bulk. Distinguishing among competing states 
by comparing their often similar condensation free ener- 
gies, requires computing these individual condensation 
energies to very high precision. This task is numeri- 
cally difficult. Here the advantages of the expression 
from Ref. which we use (see Eq. O) are obvious. 
Only the energies are explicitly needed in the sums and 
integrals on the right side, the influence of the eigen- 
functions being only indirectly felt through the relatively 
smooth function A(z). The required quantities are ob- 
tained with sufficient precision, as we shall demonstrate, 
from our numerical calculations. Some of the different 
but equivalent expressions found in the literature for the 
condensation energies or free energies^i*^, contain ex- 
plicit sums over eigenfunctions, which can lead to cumu- 
lative errors. Also, the approximate formula used for the 
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condensation energy of trilayers in Ref. 0i consisting in 
replacing, in Eq. 1)170. the pair amplitude with its spatial 
average, while plausible, is difficult to validate in general, 
and becomes inaccurate for systems involving a larger 
number of thinner layers, particularly with S widths of 
order of ho- 
using this procedure, the task then becomes feasible 
but not trivial as it involves a sum over more than 10 6 
terms. Also, to obtain the condensation energy one still 
has to subtract from the superconducting T the corre- 
sponding normal quantity T^[. It is from subtracting 
these larger quantities that the much smaller condensa- 
tion free energy is obtained. The normal state energy 
spectra used to evaluate Ttf are computed by taking 
D = in Eq. and diagonalizing the resulting ma- 
trix. The presence of interfacial scattering and mismatch 
in the Fermi wavevectors still introduces off-diagonal ma- 
trix elements. In performing the subtraction, care must 
be taken (as in fact is also the case in the bulk analytic 
case) to include exactly the same number of states in 
both calculations, rather than loosely using the same en- 
ergy cutoff: no states overall are created or lost by the 
phase transition. With all of these precautions, results 
of sufficient precision and smoothness are obtained. Re- 
sults obtained with alternative equivalent^ formulas are 
consistent but noisier. We have verified that in the lim- 
iting case of large superconducting slabs our procedure 
reproduces the analytic bulk results. 



B. Pair amplitude 

We begin by presenting some results for the pair ampli- 
tude F(z), showing how this quantity varies as a function 
of the interface scattering parameter Hb and the Fermi 
wavevector mismatch A. This is best done by means of 
three dimensional plots. In the first of these, Fig. [21 we 
show the pair amplitude (normalized to Ao) for a three 
layer SFS system, with Dp = 10 and D$ = 100, as a 
function of position and of mismatch parameter A, at 
Hb = 0. The top panel shows the results of attempting 
to find a solution of the "0" type by starting the iteration 
process with an initial guess of that form. Clearly such 
an attempt fails at small mismatch (A > 0.7) indicating 
that a solution of this type is then unstable. At larger 
mismatch, a state solution is found. Panel (b) shows 
solutions obtained by starting with an initial guess of the 
'V type: a self consistent solution of this type is then 
always obtained (note that the A scales run in different 
directions in panels (a) and (b)). For A > 0.7, the solu- 
tions of course turn out to be the same as found in panel 
(a) in the same range. Thus, we see that for small mis- 
match there is only one self consistent solution, which is 
of the 7r type, while when the mismatch is large there 
are two competing solutions and their relative stability 
becomes an issue. 

We now turn to seven layer SFSFSFS structures. In 
classifying the different possible configurations, it is con- 




FIG. 2: (Color online). The pair amplitude F(Z), normalized 
to Ao, for a three layer SFS structure, plotted as a function 
of Z = kpsz and of the mismatch parameter A, at Hb = 0. 
The direction of the A scale is different in the top and bottom 
panels. The Z = point is at the center of the structure. We 
have Ds = kFsds ~ 100 and Df = ferfp = 10. Panel (a) 
corresponds to self consistent results obtained with an initial 
guess where the junction is of the "0" type and panel (b) with 
a 'V type. In the latter case, the solution found is always 
of the 7r type, but in the former a solution of the type is 
obtained only for large mismatch (small A). We have I — 0.2 
and T = 0.01T,? here. See text for discussion. 



venient to establish a notation that envisions the seven 
layer geometry as consisting of three adjacent SFS junc- 
tions. Thus, up to a trivial reversal, we can then denote 
as "000" the structure when adjacent S layers always 
have the same sign of A(z), and as "7T7r7r" the structure 
where successive S layers alternate in sign. There are also 
two other distinct symmetric states: one in which A(z) 
has the same sign in the first two S layers, and in the last 
two it has the opposite sign, (this is labeled as the "07r0" 
configuration), and the other corresponding to the two 
outer S layers having the same sign for A(z), opposite 
to that of the two inner S layers: these are referred to 
as "7r07r" structures in this notation. We will focus our 
study on these symmetric configurations. Asymmetric 
configurations corresponding in our notation to the 7r00, 
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FIG. 3: (Color online). The normalized pair amplitude F(Z) 
for a seven layer SFSFSFS structure, plotted as in Fig. [5] 
for the same parameter values. In panels (a) and (b), the 
thickness of all S layers is the same, while in panels (c) and (d) 
the thickness of the two inner S layers is doubled to 2Ds = 200 
(see Fig.0. Panels (a) and (c) correspond to an initial guess 
of the "000" type and panels (b) and (d) to a "-kitiv" type, (see 
text). The configuration of the plotted self consistent results 
can be "000", "nnn" or "nOn" as explained in the text. 



and 7T7rO states are not forbidden, but occur very rarely 
and will be addressed only as need may arise. In Fig. El 
we repeat the plots in Fig. for seven layer structures. 
We include the cases in which all S layers are of the same 
thickness (top two panels) and the case where the thick- 
ness of the two inner S layers is doubled (bottom panels, 
see Fig. QJ. In panels (a) and (c), the initial guess is of 
the 000 type, while in (b) and (d) it is of the 7T7T7r type. 
In the case of identical S layer widths, we see that a 000 
guess (Fig. Efa)) yields a self consistent state of the same 
000 form only for larger mismatch, A < 0.5, while for 
smaller mismatch the configuration obtained is clearly of 
the 7r07r form. Thus there is a value of A where two self 
consistent solutions cross over. However, a 7T7T7r guess 
(Fig- Elk)) results in a self consistent 7T7T7r configuration 
for the whole A range. Thus, there is a clear competition 
between at least these three observed states, resulting 
from multiple minima of the free energy. Solutions of 
the 07r0 type are not displayed in this Figure but they 
will be discussed below. In the two bottom panels we see 
the same effects when the thickness of the inner S layers 
is doubled. As explained above, this describes a more 
balanced situation, since the inner layers have magnetic 
neighbors on both sides. It is evident from the figure that 
the pairing correlations are increased in the S layers. In 




FIG. 4: (Color online). The pair amplitude F(Z) for a seven 
layer structure, plotted for the same parameter values and 
conventions as in the upper panels of Fig. [5] but as a function 
of the dimensionless barrier height Hb (at A = 1), rather 
than of the mismatch. 

Fig. EJc)j there is also a noticeable shift in the crossover 
point separating the 000 and 7r07r self consistent states, 
occurring now at smaller mismatch A w 0.7. Again we 
find the competition between the various states extends 
through the entire range of A considered. 

Next we consider in Fig. El the influence of the bar- 
rier height as represented by the dimensionless parame- 
ter Hb- This figure is completely analogous to Fig.EJa) 
and (b), except for the substitution of Hb for A, which is 
set to unity (no mismatch). We find, by examining and 
comparing the panels that two solutions of the 7T7T7r and 
7r07r type exist for small barrier heights (Hb ^ 0.5), but 
when Hb becomes larger the 7T7r7r and 000 states then 
compete. Another trend which one can clearly discern is 
that the absolute values of F(Z) in the middle of the S 
layers increase with Hb- This makes sense, as at larger 
barrier heights the layers become more isolated from each 
other, and the proximity effects must in general weaken. 

To conclude this discussion of F(z) we display in Fig. El 
the pair amplitude for the same three layer system as in 
Fig. El as a function of position and of the magnetic ex- 
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FIG. 5: (Color online). The pair amplitude F(Z) for a three 
layer structure, plotted for the same parameter values and 
conventions as in Fig. [5] but as a function of the dimensionless 
magnetic polarization /, at A = 1 and Hb ~ 0. 



change parameter / in its entire range, without a barrier 
or mismatch. Careful examination of the two panels re- 
veals a rather intricate situation: a solution of the it type 
exists nearly in the entire I range (see Fig. EJb)) , the ex- 
ception being at very small /, where the magnetism be- 
comes weaker and, as one would expect, only the state 
solution is found. This requires small values of/, I < 0.1 
however. One can see that in a small neighborhood of 
/ w 0.1, as the state transitions to the ir state, the 
pair amplitude is small throughout the structure. This 
may correspond to a marked dip in the transition tem- 
perature. At larger values of I, the two types of solutions 
coexist, but there is a range around I — 0.2 where clearly 
there is only one. 

These results, which include for brevity only a very 
small subset of those obtained, are sufficient, we believe, 
to persuade the reader that although for some parame- 
ter values a unique self consistent solution exists, this is 
comparatively rare, and that in general several solutions 
of differing symmetry types can be found. These self con- 
sistent solutions correspond to local free energy minima: 
they are at least metastable. Furthermore, it is clear that 
the uniqueness or multiplicity of solutions depends in a 
complicated way not only on the geometry, but also on 
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FIG. 6: (Color online). The normalized condensation free 
energy AEo (see text) of a three layer SFS structure, plotted 
as a function of barrier height (top) and mismatch parameter 
A (bottom) for self consistent states of both the and tt types, 
as indicated. All other material parameters, geometrical val- 
ues, and temperature, are as in Fig. H The vertical arrow 
marks the end, as either Hb increases (top) or A increases 
(bottom), of the region of stability of the state in this case. 



the specific material parameter values. 



C. Condensation Free Energy: Stability 

One must, in view of the results in the previous sub- 
section, find a way to determine in each case the relative 
stability of each configuration and the global free energy 
minimum. This is achieved by computing the free energy 
of the several self consistent states, using the accurate 
numerical procedures explained earlier in this Section. 
Results for this quantity, which at the low temperature 
studied is essentially the same as the condensation en- 
ergy, are given in the next figures below. The quantity 
plotted in these figures, which we call the normalized 
A-Eq, is the condensation free energy (as calculated from 
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FIG. 7: (Color online). The normalized AEo for a seven layer 
SFSFSFS structure plotted as a function of barrier height 
(top) and mismatch parameter A (bottom) for self consistent 
states of the types indicated (see text for explanation). Ma- 
terial parameter and geometrical values are as in previous 
Figures, and all S layers are of the same thickness. The verti- 
cal arrows mark the end, as Hb decreases (top) or A increases 
(bottom, see inset) of the region of stability of a certain state 
(see text). 



Eq. (|15[) after normal state subtraction) normalized to 
N(0)Aq, which is twice the zero temperature bulk value 
(see Eg. I17I . Therefore, at the low temperatures stud- 
ied, the bulk uniform value of the quantity plotted is very 
close to -(1/2). 

In Figure we plot AEq , defined and normalized as 
explained, for a three layer SFS system. As in previous 
Figures, we have D s = 100, D F = 10 and I = 0.2. Re- 
sults for self consistent states of both the and tt type 
are plotted as indicated. The top panel shows AEo as a 
function of the barrier thickness parameter Hb at A = 1. 
The bottom panel plots the same quantity as a function 
of mismatch A at zero barrier and should be viewed in 
conjunction with Fig. [3 Looking first at the top panel, 
one sees that the zero state is stable (has nonzero con- 
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FIG. 8: (Color online). Normalized AEo for a seven layer 
SFSFSFS structure, as in Fig.Q but for the case where the 
thickness of the two inner S layers is doubled. 



densation energy) only for Hb greater than about 0.31. 
An attempt to find a solution of the type for Hb just 
below its "critical" value by using a solution of that type 
previously found for a slightly higher Hb as the starting 
guess, and iterating the self consistent process, leads after 
many iterations to a solution of the tt type. This is indi- 
cated by the vertical arrow. At larger barrier heights, 
the two states become degenerate. This makes sense 
physically: as the barriers become higher the proximity 
effect becomes less important, and the S layers behave 
more as independent superconducting slabs. The rela- 
tive phase is then immaterial. For even larger Hb we 
expect, from Eq. 1171 and the geometry, the normalized 
quantity plotted to trend, from above, toward a limit 
~ —0.5(1 — Dp/ Ds) = —0.45 and this is seen in the top 
panel. One can also see that in the region of interest 
(barriers not too high), the absolute value of the conden- 
sation energy is substantially below that of the bulk. In 
the bottom panel similar trends can be seen: in the ab- 
sence of mismatch (A = 1) only the tt state is found, and 
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FIG. 9: (Color online). The normalized AEo for a three layer 
SFS system, as a function of the thickness ds of the S layers 
(given in units of £o) at fixed Df = 10, I = 0.2 and without a 
barrier or mismatch. Results for and tt self consistent states 
are given as indicated. 



its condensation energy exhibits a somewhat oscillatory 
behavior as A decreases from unity. The state does not 
appear until A is about 0.7 and attempts (by the proce- 
dure just described) to find it lead to a 7r solution upon 
iteration (arrow). This is in agreement with the results in 
Fig. |21 For large mismatch the absolute value of AEq in- 
creases, as the S slabs become more weakly coupled, with 
a trend toward the limiting value just discussed. A very 
important difference between the top and bottom pan- 
els, however, is the crossing of the curves near A = 0.33 
in Fig. EJb). This is in effect a first order phase transi- 
tion between the tt and configurations as the mismatch 
changes. 

The results of performing the same study for a seven 
layer system with four S layers can be seen in the next 
two Figures. Figure [7| corresponds to the case where all 
the S layers have the same thickness (Ds = 100, and all 
other parameters also as in the previous figures), while 
in Fig. [SJ the thickness of the two inner S layers is dou- 
bled. Results for the four possible symmetric junction 
configurations mentioned in conjunction with Fig. [JJ are 
given, as indicated in the legends of Figs. [7| and [SJ Three 
of those configurations, 000, 7T7r7r and 7r07r have appeared 
among the results in Figs. and 0] The other configura- 
tion corresponds to the 07r0 sequence. We see that there 
are some striking differences between these examples and 
the three layer system. While in the latter case a configu- 
ration ceases to exist only when its condensation energy 
tends to zero, now configurations can become unstable 
even when, for nearby values of the relevant parameter, 
they still have a negative condensation energy. As this 
occurs, the vertical arrows in each panel indicate (an in- 
set is needed in one case for clarity) how the states trans- 
form into each other as one varies the parameters from 
the unstable to the stable region. Regardless of whether 



V 

V 

- • V T 
v v 

• 


• 

T n 


sum 


V" '" — 
• V 

• 

1 

1 1 1 1 1 1 1 1 1 


• o 

▼ n 

' i 1 ii 




- • V 


• 





0.0 0.2 0.4 06 0.& 1.D 



I 

FIG. 10: (Color online). The normalized AEo for a three 
layer SFS junction, as a function of the parameter I, for 
three different S thicknesses (as labeled) and fixed Df = 10. 
Results for and tt self consistent states are given as indicated. 



the inner layers are doubled or not, the tendency is for 
the innermost junction to remain of the same type, while 
the two outer junctions flip. Comparing Fig.[7|and Fig.[SJ 
we see that the doubling of the inner layers has a clear 
quantitative effect without having any strong qualitative 
influence. An important difference between the two cases 
is that in the first (all S layer widths equal) the two possi- 
ble states (7r07r and 7r7T7r) at zero barrier and no mismatch 
are nearly degenerate, while in the other case, the 7T7T7r 
configuration is favored. In the first case, the degener- 
acy is removed as the barrier height begins to increase, 
but A has a small effect in relative stability. In Fig. [HJ 
bottom panel, the oscillatory effect of A near the no- 
mismatch limit is seen, as in the three layer case. For 
large mismatch or barrier, the results become again de- 
generate and trend towards the appropriate limit. We 
expect these seven layer results to be at least qualita- 
tively representative of what occurs for larger values of 
Nt,: thus, states of the types 000 • • ■ 000, tttt ■ ■ ■ tttt, and 
7r00 ■ • • 007T (outer junctions one way and inner ones the 
other) should predominate for large Nl. 

It is at least of equal interest to study how the stability 
depends on the geometry. We discuss this question in the 
next four figures. First, in Fig. [5] we present results for 



13 




FIG. 11: (Color online). The normalized AEo for a three 
layer SFS system, as a function of dp, (rather than of I as in 
Fig. 1 1 ( H for four different S thicknesses (as labeled) and fixed 
7 = 0.2, H B =0, A = 1. 
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the condensation free energy of a three layer system as 
a function of ds/£,o at fixed Dp = 10, I = 0.2, Hb = 0, 
A = 1. We see that ds must be at least half a corre- 
lation length for superconductivity to be possible at all 
in this system. Convergence near that value is rather 
slow, requiring approximately 200 iterations. The super- 
conducting state then begins occurring, for this value of 
Dp, in a tt configuration only. When ds reaches £o > the tt 
state condensation energy reaches already an appreciable 
value that is consistent with that seen in the appropriate 
limits of the panels in Fig. [fJl The state is still not at- 
tainable (again, consistent with Fig. |SJ) until ds reaches 
a somewhat larger value. The condensation energies of 
the two states converge slowly toward each other upon 
increasing ds, but remain clearly non-degenerate well be- 
yond the range plotted. The small breaks in the state 
curve correspond to specific S widths that permit only 
the 7T state. 

The behavior seen in Fig. |5| depends strongly on /. 
This dependence is displayed in Fig. ^| where we show 
the normalized AEq for the same system, as a function of 
I, for three different values of ds- For the value I = 0.2, 
the results shown are consistent with Fig. ^ including 
the nonexistence of the state at ds = £o- We now 
see, however, that it is not always the tt state which is 
favored, but that the difference in condensation energies 
is an oscillatory function of I. This of course reflects 
that whether the or the tt state is preferred depends, all 
other things being equal, on the relation between Dp and 
(fcj — fc|) _1 , and this quantity (see Eq. ijTol) ) depends on 
I. At intermediate values of I (centered around I = 0.5) 
the zero state is favored, and when I becomes very small 
the tt state ceases to exist altogether (this can be seen also 
by examination of Fig. [5J. As I — > the condensation 
energy of the state remains somewhat above the bulk 



FIG. 12: (Color online). The normalized condensation free 
energy for a seven layer SFSFSFS system, as a function 
of D F , for fixed I = 0.2, d s = £o, A = 1 and H B = 0.5. 
Results for the four possible symmetric self consistent states 
are given, as indicated, for both the cases where all S layers 
are identical (top) and where the thickness of the inner ones 
is doubled (bottom). Lines are guides to the eye. Breaks (top 
panel) indicate regions where a certain configuration is not 
found. 



value and, as one would expect, decreases slightly with 
increasing ds- At larger values of /, the absolute values of 
AEq increase with ds and on the average decrease slowly 
with /. 

The oscillations in Fig. as a function of I at con- 
stant dp can also be displayed by considering results as a 
function of dp at constant I. We do so in Fig. 1111 where 
we plot, for a three layer SFS system, AEo as a function 
of Dp at constant I = 0.2 for four values of ds/£,o- One 
sees again that for this value of I the state does not 
exist at ds = £o & n d Dp = 10 but that it appears at 
larger values of ds/Co- The damped oscillatory behavior 
is quite evident. At larger values of dp the condensa- 
tion energies of the two states trend towards a common 
value that increases in absolute value with ds- At a very 
small value of dp, which depends on ds, the tt state be- 
gins to vanish, and the condensation free energy of the 
state tends then towards the bulk value. All of this is 
consistent with simple physical arguments. 

In Fig.^Jwe extend the results of Fig.^Jto the seven 
layer system. In this case we consider only one value of 
ds (ds = £o) but include a finite barrier thickness, Hb = 
0.5. The finite barrier allows for the possibility of more 
distinct states coexisting (see Figs. 17181 . We consider 
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both the cases where all S layers thicknesses are equal 
(top panel) and the case where the inner ones are doubled 
(bottom panel). All possible symmetric self consistent 
states were studied, as indicated in the figure. In contrast 
with the three layer example with no barrier, in the seven 
layer cases with Hb = 0.5 all of the four symmetric states 
(000, 7T7r7r, ttOtt, OttO) are at least metastable over a range 
of dp, even at ds = £o- ln the top panel we see however 
that only the 07r0 state is stable over the whole dp range. 
The 7r7T7r state reverts to the 07r0 state in the range 1.6 < 
dp < 4.2, while the 7r07r state reverts to 000 state for 
1.8 < dp < 2.4. The 000 state is unstable for much of 
the range for 6 < dp < 12. It appears that in this range 
the 000 state is sufficiently close to a crossover (see e.g. 
Fig. 01 that attempts to find it sometimes converge to 
an asymmetric 7r00 state, rather than the expected ttOtt. 
In these cases the number of iterations to convergence is 
substantially increased, as the order parameter attempts 
to readjust its profile. For the situation where the inner 
S layers are twice the width of the outer ones, we see 
(bottom panel) that all four symmetric configurations are 
either stable or metastable for the whole dp range. This 
is consistent with Fig.[7| where at Hb = 0.5 all four states 
are present simultaneously. The condensation energy is 
of course lower than in the previous case, due to the 
increased pairing correlations associated with the thicker 
S slabs. For both geometries oscillations arising from the 
scattering potential lead to deviations from the estimated 
periodicity determined by (fcf — fc|) _1 . For sufficiently 
large dp the difference in energies becomes small. One 
can infer from these results than in superlattices with 
realistic oxide barriers, where as the number of layers 
increases a larger number of nontrivial possible states 
arise, the number of local free energy minima that can 
coexist will increase. 



D. DOS 

We conclude this section by presenting some results 
for the DOS, as it can be experimentally measured. The 
results given are in all cases for the quantity N a (z, e) de- 
fined in Eq. (|14f) . summed over cr, and integrated over 
a distance of one coherence length from the edge of the 
sample. We normalize our results to the corresponding 
value for the normal state of the superconducting mate- 
rial, and the energies to the bulk value of the gap, A . 
All parameters not otherwise mentioned are set to their 
"default" values outlined at the beginning of this section. 

In the top panel of Fig. ^| we present results for an 
SFS trilayer, for two contrasting values of the mismatch 
parameter A. As usual, results labeled as "0" and 'V 
are for the case where the self consistent states plot- 
ted are of these types. The and tt state curves cor- 
responding to A = 0.67, where (see Fig. EJ) the state is 
barely metastable, have clearly distinct signatures, with 
a smaller gap opening for the state, and consequently 
more subgap quasiparticle states. Therefore, as can be 




FIG. 13: (Color online). The normalized DOS (see text) 
for a SFS trilayer, plotted as a function of the energy (in 
units of the bulk zero temperature gap). Results for both 
and n self consistent states are given, as indicated. In the 
top panel, the DOS is shown at Hb = for two different 
values of the mismatch parameter, A = 0.2, and 0.67, the 
latter being a case for which the state is nearly unstable 
(see Fig.|5J. The bottom panel shows the DOS profile in the 
absence of mismatch (A = 1), but with the interface scattering 
parameter Hb taking on the two values shown, chosen on 
similar criteria as the A values (see text). 



seen in Fig. [3J when there is little mismatch the pair am- 
plitude is relatively large in the F layer. When there 
is more mismatch the proximity effect weakens, the gap 
opens and large peaks form which progressively become 
more BCS-like as the mismatch increases (A decreases) to 
A = 0.2. This progression takes a different form for the 
and 7r states, as one can see by careful comparison of the 
curves. In the bottom panel we demonstrate the effect of 
the barrier height. The results are displayed as in the top 
panel, but with the dimensionless height Hb taking the 
place of the mismatch parameter. This Figure should be 
viewed in conjunction with the top panel of Fig. [SJ One 
of the values of Hb chosen (Hb = 0.32) is again such that 
the state is barely metastable, while for the other value 
the and 7r states have similar condensation energies. 
For Hb close to 0.31 there is a marked contrast between 
the two plots, with the gap clearly opening wider, and 
containing more structure for the tt state. At larger Hb 
the gap becomes larger in both cases, with the BCS-like 
peaks increasing in height. Thus, as the barriers becomes 
larger, one is dealing with nearly independent supercon- 
ducting slabs and the plots become more similar. The 
largest difference therefore, occurs, as for the mismatch, 
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FIG. 14: (Color online). The normalized local DOS for a 
seven layer system versus the dimensionless energy and mis- 
match parameter A. The panels are arranged as those in 
Fig- 01 ( a ) an d (b) show results obtained when all S layers 
have the same thickness, while in (c) and (d) that thickness 
is doubled. Panels (a) and (c) correspond to starting the it- 
eration procedure with a 000 order parameter, while panels 
(b) and (d) are obtained by starting with a -Kirn junction, in 
which case the self-consistent solution is also always of this 
type. 



at the intermediate values more likely to be found exper- 
imentally. 

It is also illustrative to display the DOS in a three 
dimensional format. Thus, in Fig. 1141 the normalized 
DOS is presented as a function of e/A and A for the 
seven layer structures previously discussed. This permits 
a much more extensive range of A values to be exam- 
ined. The panel arrangement is the same as in Fig. El in 
the top two panels, all S layers have the same thickness 
(ds — £o) while in the bottom ones the thickness of the 
inner ones is doubled. On the right column (panels (b) 
and (d)) the self-consistent results shown were obtained 
by starting the iteration process with an initial guess of 
the 7T7r7r type, in which case, as one can see from Fig. UJ 
and the bottom panels of Figs. and [SJ the solution is 
also of the 7T7r7r type. The left column panels were ob- 
tained from an initial guess of the 000 type and therefore 
correspond, as one can see from Figs. [3J[71 and [SJ to solu- 
tions of the 000 type for A < 0.5 (panel (a)) and A < 0.7 
(panel (c)), and to solutions of the 07r0 type otherwise. 
There is no overlap among the results in the four panels. 
One sees again that as mismatch increases the BCS like 
peaks become more prominent and move out, trending 
towards their bulk positions, while the gap opens. Fur- 
thermore, the 7T7T7r results for the doubled ds inner layers 
are smoother and show clearly a broader gap throughout, 
due to the extended S geometry. The DOS behavior ob- 



served here is entirely consistent with the condensation 
energy results found. 



IV. CONCLUSIONS 

In summary, we have found self consistent solutions to 
the microscopic BdG equations for SFS and SFSFSFS 
structures, for a wide range of parameter values. We 
have shown that, in most cases, several such self consis- 
tent solutions coexist, with differing spatial dependence 
of the pair potential A(r) and the pair amplitude F(r). 
Thus, there can be in general competing local minima of 
the free energy. Determining their relative stability re- 
quires the computation of their respective condensation 
free energies, which we have done by using an efficient, 
accurate approach that does not involve the quasiparticle 
amplitudes directly, and requires only the eigenenergies 
and the pair potential. 

For SFS trilayers (single junctions), we found that 
both 7r and junction states exist for a range of values 
of the relevant parameters. We have displayed results 
for the pair amplitude, which give insight into the su- 
perconducting correlations, and for the condensation free 
energies of each configuration, to determine the true equi- 
librium state. We have shown that a transition (which 
is in effect of first order in parameter space) can occur 
between the and 7r states for a critical value of A. The 
difference in condensation energies between the two pos- 
sible states exhibits oscillations as a function of I. This 
behavior is strongly dependent on the width of the S lay- 
ers. For ds equal to one coherence length £oj there exists 
a range of I in which either a or tt state survives, but 
not both. Increasing the S width by about £o/2 restores 
the coexistence of both states. 

Several interesting phenomena arise when one explores 
the geometrical parameters of trilayer structures. For a 
fixed ferromagnet width dp, and parameters values that 
lead to a 7r state, we found (see Fig. [5J that the 7r con- 
figuration remains the ground state of the system as ds 
varies. The 7r state first appears at small ds (ds ~ Co/2), 
and then its condensation energy declines monotonically 
towards the bulk limit. The metastable state begins at 
larger ds ~ Co, and its condensation energy declines also 
slowly over the range of ds studied. The other relevant 
length that was considered is dp- The condensation en- 
ergy is, for both states, an oscillatory function of dp. The 
oscillations become better defined, and the possibility of 
both and 7r states coexisting increases, at larger ds- 
As expected, we find that the condensation energy has 
very similar properties when either dp or / varies. The 
period approximately agrees with the estimate given by 
(fef — fc^) _1 , which governs the oscillations of the pair am- 
plitude and in general, of other single particle quantities. 
The results for the DOS reflect the crossover from a state 
with populated subgap peaks to a nearly gapped BCS- 
like behavior as A is decreased or the barrier height Hb 
increased. Signatures that may help to experimentally 
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identify the and 7r configurations were seen. 

As the number of layers increases, so does the number 
of competing stable and metastable junction configura- 
tions. We considered two types of seven layers structures, 
and found that doubling the width of the inner S lay- 
ers (which are bounded on each side by ferromagnets), 
resulted typically in different quasiparticle spectra and 
pair amplitudes, compared to the situation when all S 
layers have the same ds- For large mismatch or barrier 
strength, the phase of the pair amplitude in each layer 
is independent, and configurations are nearly degener- 
ate, but as each of these parameters diminishes there 
is a crossing over to a situation where the free energies 
of each configuration are well-separated. At certain val- 
ues of A and Hb, some configurations become unstable. 
These values are different depending on the type of sys- 
tem (single or double inner layers). For fixed / and dp 
we found that if a state is stable at no mismatch and 
zero barrier, then it remains at least metastable over a 
very wide range of A and Hb values. Our results showed 
that self-consistency cannot be neglected as the number 
of layers increases, due to the nontrivial and intricate 
spatial variations in F(r) that become possible. 

For seven layers, we studied in detail the condensa- 
tion free energies of the four symmetric junction states, 
000, 7T7r7r, 7r07r and 07r0, in the previously introduced no- 
tation. We first investigated the stable states as a func- 
tion of A and Hb ■ In contrast to the three layer system, 
we found that states could become unstable even when 
the condensation energy did not tend to zero for nearby 
values of the relevant parameters. For double width in- 
ner layered structures, we found a greater spread in the 
free energies between the four states, and the instabil- 
ity found in certain cases for the 000 and 07r0 states was 
shifted in A and Hb, in agreement with the pair ampli- 
tude results. It is reasonable to assume that these results 
are representative of what occurs for superlattices. We 
again found transitions upon varying A and the number 
and sequence of the transitions is now more intricate (see 
Figs. and IHJ). The analysis of the geometrical properties 
revealed that scattering at the interfaces modifies the ex- 
pected damped oscillatory behavior of the condensation 
energy as a function of dp- In effect, the barriers intro- 
duce significant atomic scale oscillations that smear the 



periodicity. This underlines the importance of a micro- 
scopic approach for the investigation of nanostrucutres. 
As with the A dependence, we also showed that the global 
minima in the free energy is different for the two struc- 
tures as dp changes. The configuration of the ground 
state of the system with S layers of uniform width was 
more variable in parameter space, compared to when the 
inner layers are doubled, (compare Figs. [7] and |8J). Fi- 
nally, we calculated the DOS, to illustrate and compare 
the differences in the spectra for the two different seven 
layer geometries. Of the two, the energy gap for the single 
inner layer case, and 7T7r7r stable state, contains more sub- 
gap states, for the specific examples plotted fFigs. Upland 
1141 due to stronger pair breaking effects of the exchange 
field. These states fill in increasingly with decreased mis- 
match. 

Our results were obtained in the clean limit, which is 
appropriate for the relatively thin structures envisioned 
here. Furthermore, as shown in Ref. 28 in conjunction 
with realistic comparison with experiments^ the influ- 
ence of impurities can be taken into account by replac- 
ing the clean value of £o with an effective one. A sep- 
arate important issue is that of the free energy barri- 
ers separating the different free energy minima we have 
found, and hence to which degree are metastable states 
long lived. Our method cannot directly answer this ques- 
tion, but from the macroscopic symmetry differences in 
the pair amplitude structure of the different states one 
would have to conclude that the barriers are high and 
the metastable states could be very long lived. We ex- 
pect that the transitions found here in parameter space 
at constant temperature will be reflected in actual first 
order phase transitions as a function of T. Such transi- 
tions would presumably be very hysteretic. We hope to 
examine this question in the future. 
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